home *** CD-ROM | disk | FTP | other *** search
/ The Datafile PD-CD 1 Issue 2 / PDCD-1 - Issue 02.iso / _utilities / utilities / 001 / meschach / !Meschach / c / ztorture < prev    next >
Text File  |  1994-08-04  |  20KB  |  722 lines
  1.  
  2. /**************************************************************************
  3. **
  4. ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  5. **
  6. **                 Meschach Library
  7. ** 
  8. ** This Meschach Library is provided "as is" without any express 
  9. ** or implied warranty of any kind with respect to this software. 
  10. ** In particular the authors shall not be liable for any direct, 
  11. ** indirect, special, incidental or consequential damages arising 
  12. ** in any way from use of the software.
  13. ** 
  14. ** Everyone is granted permission to copy, modify and redistribute this
  15. ** Meschach Library, provided:
  16. **  1.  All copies contain this copyright notice.
  17. **  2.  All modified copies shall carry a notice stating who
  18. **      made the last modification and the date of such modification.
  19. **  3.  No charge is made for this software or works derived from it.  
  20. **      This clause shall not be construed as constraining other software
  21. **      distributed on the same medium as this software, nor is a
  22. **      distribution fee considered a charge.
  23. **
  24. ***************************************************************************/
  25.  
  26.  
  27. /*
  28.     This file contains a series of tests for the Meschach matrix
  29.     library, complex routines
  30. */
  31.  
  32. static char rcsid[] = "$Id: $";
  33.  
  34. #include    <stdio.h>
  35. #include    <math.h>
  36. #include     "zmatrix2.h"
  37. #include        "matlab.h"
  38.  
  39.  
  40. #define    errmesg(mesg)    printf("Error: %s error: line %d\n",mesg,__LINE__)
  41. #define notice(mesg)    printf("# Testing %s...\n",mesg);
  42.  
  43. /* extern    int    malloc_chain_check(); */
  44. /* #define MEMCHK() if ( malloc_chain_check(0) ) \
  45. { printf("Error in malloc chain: \"%s\", line %d\n", \
  46.      __FILE__, __LINE__); exit(0); } */
  47. #define    MEMCHK() 
  48.  
  49. /* cmp_perm -- returns 1 if pi1 == pi2, 0 otherwise */
  50. int    cmp_perm(pi1, pi2)
  51. PERM    *pi1, *pi2;
  52. {
  53.     int        i;
  54.  
  55.     if ( ! pi1 || ! pi2 )
  56.     error(E_NULL,"cmp_perm");
  57.     if ( pi1->size != pi2->size )
  58.     return 0;
  59.     for ( i = 0; i < pi1->size; i++ )
  60.     if ( pi1->pe[i] != pi2->pe[i] )
  61.         return 0;
  62.     return 1;
  63. }
  64.  
  65. /* px_rand -- generates sort-of random permutation */
  66. PERM    *px_rand(pi)
  67. PERM    *pi;
  68. {
  69.     int        i, j, k;
  70.  
  71.     if ( ! pi )
  72.     error(E_NULL,"px_rand");
  73.  
  74.     for ( i = 0; i < 3*pi->size; i++ )
  75.     {
  76.     j = (rand() >> 8) % pi->size;
  77.     k = (rand() >> 8) % pi->size;
  78.     px_transp(pi,j,k);
  79.     }
  80.  
  81.     return pi;
  82. }
  83. #ifdef RISC_OS
  84. #define SAVE_FILE "asx521_mat"
  85. #else
  86. #define    SAVE_FILE    "asx5213a.mat"
  87. #endif
  88. #define    MATLAB_NAME    "alpha"
  89. char    name[81] = MATLAB_NAME;
  90.  
  91. void    main(argc, argv)
  92. int    argc;
  93. char    *argv[];
  94. {
  95.     ZVEC     *x = ZVNULL, *y = ZVNULL, *z = ZVNULL, *u = ZVNULL;
  96.     ZVEC    *diag = ZVNULL;
  97.     PERM    *pi1 = PNULL, *pi2 = PNULL, *pivot = PNULL;
  98.     ZMAT    *A = ZMNULL, *B = ZMNULL, *C = ZMNULL, *D = ZMNULL,
  99.     *Q = ZMNULL;
  100.     complex    ONE;
  101.     complex    z1, z2, z3;
  102.     Real    cond_est, s1, s2, s3;
  103.     int        i, seed;
  104.     FILE    *fp;
  105.     char    *cp;
  106.  
  107.  
  108.     mem_info_on(TRUE);
  109.  
  110.     setbuf(stdout,(char *)NULL);
  111.  
  112.     seed = 1111;
  113.     if ( argc > 2 )
  114.     {
  115.     printf("usage: %s [seed]\n",argv[0]);
  116.     exit(0);
  117.     }
  118.     else if ( argc == 2 )
  119.     sscanf(argv[1], "%d", &seed);
  120.  
  121.     /* set seed for rand() */
  122.     smrand(seed);
  123.  
  124.     /* print out version information */
  125.     m_version();
  126.  
  127.     printf("# Meschach Complex numbers & vectors torture test\n\n");
  128.     printf("# grep \"^Error\" the output for a listing of errors\n");
  129.     printf("# Don't panic if you see \"Error\" appearing; \n");
  130.     printf("# Also check the reported size of error\n");
  131.     printf("# This program uses randomly generated problems and therefore\n");
  132.     printf("# may occasionally produce ill-conditioned problems\n");
  133.     printf("# Therefore check the size of the error compared with MACHEPS\n");
  134.     printf("# If the error is within 1000*MACHEPS then don't worry\n");
  135.     printf("# If you get an error of size 0.1 or larger there is \n");
  136.     printf("# probably a bug in the code or the compilation procedure\n\n");
  137.     printf("# seed = %d\n",seed);
  138.  
  139.     printf("\n");
  140.  
  141.     mem_stat_mark(1);
  142.  
  143.     notice("complex arithmetic & special functions");
  144.  
  145.     ONE = zmake(1.0,0.0);
  146.     printf("# ONE = ");    z_output(ONE);
  147.     z1.re = mrand();    z1.im = mrand();
  148.     z2.re = mrand();    z2.im = mrand();
  149.     z3 = zadd(z1,z2);
  150.     if ( fabs(z1.re+z2.re-z3.re) + fabs(z1.im+z2.im-z3.im) > 10*MACHEPS )
  151.     errmesg("zadd");
  152.     z3 = zsub(z1,z2);
  153.     if ( fabs(z1.re-z2.re-z3.re) + fabs(z1.im-z2.im-z3.im) > 10*MACHEPS )
  154.     errmesg("zadd");
  155.     z3 = zmlt(z1,z2);
  156.     if ( fabs(z1.re*z2.re - z1.im*z2.im - z3.re) +
  157.      fabs(z1.im*z2.re + z1.re*z2.im - z3.im) > 10*MACHEPS )
  158.     errmesg("zmlt");
  159.     s1 = zabs(z1);
  160.     if ( fabs(s1*s1 - (z1.re*z1.re+z1.im*z1.im)) > 10*MACHEPS )
  161.     errmesg("zabs");
  162.     if ( zabs(zsub(z1,zmlt(z2,zdiv(z1,z2)))) > 10*MACHEPS ||
  163.      zabs(zsub(ONE,zdiv(z1,zmlt(z2,zdiv(z1,z2))))) > 10*MACHEPS )
  164.     errmesg("zdiv");
  165.  
  166.     z3 = zsqrt(z1);
  167.     if ( zabs(zsub(z1,zmlt(z3,z3))) > 10*MACHEPS )
  168.     errmesg("zsqrt");
  169.     if ( zabs(zsub(z1,zlog(zexp(z1)))) > 10*MACHEPS )
  170.     errmesg("zexp/zlog");
  171.     
  172.  
  173.     printf("# Check: MACHEPS = %g\n",MACHEPS);
  174.     /* allocate, initialise, copy and resize operations */
  175.     /* ZVEC */
  176.     notice("vector initialise, copy & resize");
  177.     x = zv_get(12);
  178.     y = zv_get(15);
  179.     z = zv_get(12);
  180.     zv_rand(x);
  181.     zv_rand(y);
  182.     z = zv_copy(x,z);
  183.     if ( zv_norm2(zv_sub(x,z,z)) >= MACHEPS )
  184.     errmesg("ZVEC copy");
  185.     zv_copy(x,y);
  186.     x = zv_resize(x,10);
  187.     y = zv_resize(y,10);
  188.     if ( zv_norm2(zv_sub(x,y,z)) >= MACHEPS )
  189.     errmesg("ZVEC copy/resize");
  190.     x = zv_resize(x,15);
  191.     y = zv_resize(y,15);
  192.     if ( zv_norm2(zv_sub(x,y,z)) >= MACHEPS )
  193.     errmesg("VZEC resize");
  194.  
  195.     /* ZMAT */
  196.     notice("matrix initialise, copy & resize");
  197.     A = zm_get(8,5);
  198.     B = zm_get(3,9);
  199.     C = zm_get(8,5);
  200.     zm_rand(A);
  201.     zm_rand(B);
  202.     C = zm_copy(A,C);
  203.     if ( zm_norm_inf(zm_sub(A,C,C)) >= MACHEPS )
  204.     errmesg("ZMAT copy");
  205.     zm_copy(A,B);
  206.     A = zm_resize(A,3,5);
  207.     B = zm_resize(B,3,5);
  208.     if ( zm_norm_inf(zm_sub(A,B,C)) >= MACHEPS )
  209.     errmesg("ZMAT copy/resize");
  210.     A = zm_resize(A,10,10);
  211.     B = zm_resize(B,10,10);
  212.     if ( zm_norm_inf(zm_sub(A,B,C)) >= MACHEPS )
  213.     errmesg("ZMAT resize");
  214.  
  215.     MEMCHK();
  216.  
  217.     /* PERM */
  218.     notice("permutation initialise, inverting & permuting vectors");
  219.     pi1 = px_get(15);
  220.     pi2 = px_get(12);
  221.     px_rand(pi1);
  222.     zv_rand(x);
  223.     px_zvec(pi1,x,z);
  224.     y = zv_resize(y,x->dim);
  225.     pxinv_zvec(pi1,z,y);
  226.     if ( zv_norm2(zv_sub(x,y,z)) >= MACHEPS )
  227.     errmesg("PERMute vector");
  228.  
  229.     /* testing catch() etc */
  230.     notice("error handling routines");
  231.     catch(E_NULL,
  232.       catchall(zv_add(ZVNULL,ZVNULL,ZVNULL);
  233.              errmesg("tracecatch() failure"),
  234.              printf("# tracecatch() caught error\n");
  235.              error(E_NULL,"main"));
  236.                  errmesg("catch() failure"),
  237.       printf("# catch() caught E_NULL error\n"));
  238.  
  239.     /* testing inner products and v_mltadd() etc */
  240.     notice("inner products and linear combinations");
  241.     u = zv_get(x->dim);
  242.     zv_rand(u);
  243.     zv_rand(x);
  244.     zv_resize(y,x->dim);
  245.     zv_rand(y);
  246.     zv_mltadd(y,x,zneg(zdiv(zin_prod(x,y),zin_prod(x,x))),z);
  247.     if ( zabs(zin_prod(x,z)) >= 5*MACHEPS*x->dim )
  248.     {
  249.     errmesg("zv_mltadd()/zin_prod()");
  250.     printf("# error norm = %g\n", zabs(zin_prod(x,z)));
  251.     }
  252.  
  253.     z1 = zneg(zdiv(zin_prod(x,y),zmake(zv_norm2(x)*zv_norm2(x),0.0)));
  254.     zv_mlt(z1,x,u);
  255.     zv_add(y,u,u);
  256.     if ( zv_norm2(zv_sub(u,z,u)) >= MACHEPS*x->dim )
  257.     {
  258.     errmesg("zv_mlt()/zv_norm2()");
  259.     printf("# error norm = %g\n", zv_norm2(u));
  260.     }
  261.  
  262. #ifdef ANSI_C
  263.     zv_linlist(u,x,z1,y,ONE,VNULL);
  264.     if ( zv_norm2(zv_sub(u,z,u)) >= MACHEPS*x->dim )
  265.     errmesg("zv_linlist()");
  266. #endif
  267. #ifdef VARARGS
  268.     zv_linlist(u,x,z1,y,ONE,VNULL);
  269.     if ( zv_norm2(zv_sub(u,z,u)) >= MACHEPS*x->dim )
  270.     errmesg("zv_linlist()");
  271. #endif
  272.  
  273.     MEMCHK();
  274.  
  275.     /* vector norms */
  276.     notice("vector norms");
  277.     x = zv_resize(x,12);
  278.     zv_rand(x);
  279.     for ( i = 0; i < x->dim; i++ )
  280.     if ( zabs(v_entry(x,i)) >= 0.7 )
  281.         v_set_val(x,i,ONE);
  282.         else
  283.         v_set_val(x,i,zneg(ONE));
  284.     s1 = zv_norm1(x);
  285.     s2 = zv_norm2(x);    
  286.     s3 = zv_norm_inf(x);
  287.     if ( fabs(s1 - x->dim) >= MACHEPS*x->dim ||
  288.      fabs(s2 - sqrt((double)(x->dim))) >= MACHEPS*x->dim ||
  289.      fabs(s3 - 1.0) >= MACHEPS )
  290.     errmesg("zv_norm1/2/_inf()");
  291.  
  292.     /* test matrix multiply etc */
  293.     notice("matrix multiply and invert");
  294.     A = zm_resize(A,10,10);
  295.     B = zm_resize(B,10,10);
  296.     zm_rand(A);
  297.     zm_inverse(A,B);
  298.     zm_mlt(A,B,C);
  299.     for ( i = 0; i < C->m; i++ )
  300.     m_set_val(C,i,i,zsub(m_entry(C,i,i),ONE));
  301.     if ( zm_norm_inf(C) >= MACHEPS*zm_norm_inf(A)*zm_norm_inf(B)*5 )
  302.     errmesg("zm_inverse()/zm_mlt()");
  303.  
  304.     MEMCHK();
  305.  
  306.     /* ... and adjoints */
  307.     notice("adjoints and adjoint-multiplies");
  308.     zm_adjoint(A,A);    /* can do square matrices in situ */
  309.     zmam_mlt(A,B,C);
  310.     for ( i = 0; i < C->m; i++ )
  311.     m_set_val(C,i,i,zsub(m_entry(C,i,i),ONE));
  312.     if ( zm_norm_inf(C) >= MACHEPS*zm_norm_inf(A)*zm_norm_inf(B)*5 )
  313.     errmesg("zm_adjoint()/zmam_mlt()");
  314.     zm_adjoint(A,A);
  315.     zm_adjoint(B,B);
  316.     zmma_mlt(A,B,C);
  317.     for ( i = 0; i < C->m; i++ )
  318.     m_set_val(C,i,i,zsub(m_entry(C,i,i),ONE));
  319.     if ( zm_norm_inf(C) >= MACHEPS*zm_norm_inf(A)*zm_norm_inf(B)*5 )
  320.     errmesg("zm_adjoint()/zmma_mlt()");
  321.     zsm_mlt(zmake(3.71,2.753),B,B);
  322.     zmma_mlt(A,B,C);
  323.     for ( i = 0; i < C->m; i++ )
  324.     m_set_val(C,i,i,zsub(m_entry(C,i,i),zmake(3.71,-2.753)));
  325.     if ( zm_norm_inf(C) >= MACHEPS*zm_norm_inf(A)*zm_norm_inf(B)*5 )
  326.     errmesg("szm_mlt()/zmma_mlt()");
  327.     zm_adjoint(B,B);
  328.     zsm_mlt(zdiv(ONE,zmake(3.71,-2.753)),B,B);
  329.  
  330.     MEMCHK();
  331.  
  332.     /* ... and matrix-vector multiplies */
  333.     notice("matrix-vector multiplies");
  334.     x = zv_resize(x,A->n);
  335.     y = zv_resize(y,A->m);
  336.     z = zv_resize(z,A->m);
  337.     u = zv_resize(u,A->n);
  338.     zv_rand(x);
  339.     zv_rand(y);
  340.     zmv_mlt(A,x,z);
  341.     z1 = zin_prod(y,z);
  342.     zvm_mlt(A,y,u);
  343.     z2 = zin_prod(u,x);
  344.     if ( zabs(zsub(z1,z2)) >= (MACHEPS*x->dim)*x->dim )
  345.     {
  346.     errmesg("zmv_mlt()/zvm_mlt()");
  347.     printf("# difference between inner products is %g\n",
  348.            zabs(zsub(z1,z2)));
  349.     }
  350.     zmv_mlt(B,z,u);
  351.     if ( zv_norm2(zv_sub(u,x,u)) >= MACHEPS*zm_norm_inf(A)*zm_norm_inf(B)*5 )
  352.     errmesg("zmv_mlt()/zvm_mlt()");
  353.  
  354.     MEMCHK();
  355.  
  356.     /* get/set row/col */
  357.     notice("getting and setting rows and cols");
  358.     x = zv_resize(x,A->n);
  359.     y = zv_resize(y,B->m);
  360.     x = zget_row(A,3,x);
  361.     y = zget_col(B,3,y);
  362.     if ( zabs(zsub(_zin_prod(x,y,0,Z_NOCONJ),ONE)) >=
  363.     MACHEPS*zm_norm_inf(A)*zm_norm_inf(B)*5 )
  364.     errmesg("zget_row()/zget_col()");
  365.     zv_mlt(zmake(-1.0,0.0),x,x);
  366.     zv_mlt(zmake(-1.0,0.0),y,y);
  367.     zset_row(A,3,x);
  368.     zset_col(B,3,y);
  369.     zm_mlt(A,B,C);
  370.     for ( i = 0; i < C->m; i++ )
  371.     m_set_val(C,i,i,zsub(m_entry(C,i,i),ONE));
  372.     if ( zm_norm_inf(C) >= MACHEPS*zm_norm_inf(A)*zm_norm_inf(B)*5 )
  373.     errmesg("zset_row()/zset_col()");
  374.  
  375.     MEMCHK();
  376.  
  377.     /* matrix norms */
  378.     notice("matrix norms");
  379.     A = zm_resize(A,11,15);
  380.     zm_rand(A);
  381.     s1 = zm_norm_inf(A);
  382.     B = zm_adjoint(A,B);
  383.     s2 = zm_norm1(B);
  384.     if ( fabs(s1 - s2) >= MACHEPS*A->m )
  385.     errmesg("zm_norm1()/zm_norm_inf()");
  386.     C = zmam_mlt(A,A,C);
  387.     z1.re = z1.im = 0.0;
  388.     for ( i = 0; i < C->m && i < C->n; i++ )
  389.     z1 = zadd(z1,m_entry(C,i,i));
  390.     if ( fabs(sqrt(z1.re) - zm_norm_frob(A)) >= MACHEPS*A->m*A->n )
  391.     errmesg("zm_norm_frob");
  392.  
  393.     MEMCHK();
  394.     
  395.     /* permuting rows and columns */
  396.     /******************************
  397.     notice("permuting rows & cols");
  398.     A = zm_resize(A,11,15);
  399.     B = zm_resize(B,11,15);
  400.     pi1 = px_resize(pi1,A->m);
  401.     px_rand(pi1);
  402.     x = zv_resize(x,A->n);
  403.     y = zmv_mlt(A,x,y);
  404.     px_rows(pi1,A,B);
  405.     px_zvec(pi1,y,z);
  406.     zmv_mlt(B,x,u);
  407.     if ( zv_norm2(zv_sub(z,u,u)) >= MACHEPS*A->m )
  408.     errmesg("px_rows()");
  409.     pi1 = px_resize(pi1,A->n);
  410.     px_rand(pi1);
  411.     px_cols(pi1,A,B);
  412.     pxinv_zvec(pi1,x,z);
  413.     zmv_mlt(B,z,u);
  414.     if ( zv_norm2(zv_sub(y,u,u)) >= MACHEPS*A->n )
  415.     errmesg("px_cols()");
  416.     ******************************/
  417.  
  418.     MEMCHK();
  419.  
  420.     /* MATLAB save/load */
  421.     notice("MATLAB save/load");
  422.     A = zm_resize(A,12,11);
  423.     if ( (fp=fopen(SAVE_FILE,"w")) == (FILE *)NULL )
  424.     printf("Cannot perform MATLAB save/load test\n");
  425.     else
  426.     {
  427.     zm_rand(A);
  428.     zm_save(fp, A, name);
  429.     fclose(fp);
  430.     if ( (fp=fopen(SAVE_FILE,"r")) == (FILE *)NULL )
  431.         printf("Cannot open save file \"%s\"\n",SAVE_FILE);
  432.     else
  433.     {
  434.         ZM_FREE(B);
  435.         B = zm_load(fp,&cp);
  436.         if ( strcmp(name,cp) || zm_norm1(zm_sub(A,B,C)) >=
  437.          MACHEPS*A->m )
  438.         {
  439.         errmesg("zm_load()/zm_save()");
  440.         printf("# orig. name = %s, restored name = %s\n", name, cp);
  441.         printf("# orig. A =\n");    zm_output(A);
  442.         printf("# restored A =\n");    zm_output(B);
  443.         }
  444.     }
  445.     }
  446.  
  447.     MEMCHK();
  448.  
  449.     /* Now, onto matrix factorisations */
  450.     A = zm_resize(A,10,10);
  451.     B = zm_resize(B,A->m,A->n);
  452.     zm_copy(A,B);
  453.     x = zv_resize(x,A->n);
  454.     y = zv_resize(y,A->m);
  455.     z = zv_resize(z,A->n);
  456.     u = zv_resize(u,A->m);
  457.     zv_rand(x);
  458.     zmv_mlt(B,x,y);
  459.     z = zv_copy(x,z);
  460.  
  461.     notice("LU factor/solve");
  462.     pivot = px_get(A->m);
  463.     zLUfactor(A,pivot);
  464.     tracecatch(zLUsolve(A,pivot,y,x),"main");
  465.     tracecatch(cond_est = zLUcondest(A,pivot),"main");
  466.     printf("# cond(A) approx= %g\n", cond_est);
  467.     if ( zv_norm2(zv_sub(x,z,u)) >= MACHEPS*zv_norm2(x)*cond_est)
  468.     {
  469.     errmesg("zLUfactor()/zLUsolve()");
  470.     printf("# LU solution error = %g [cf MACHEPS = %g]\n",
  471.            zv_norm2(zv_sub(x,z,u)), MACHEPS);
  472.     }
  473.  
  474.  
  475.     zv_copy(y,x);
  476.     tracecatch(zLUsolve(A,pivot,x,x),"main");
  477.     tracecatch(cond_est = zLUcondest(A,pivot),"main");
  478.     if ( zv_norm2(zv_sub(x,z,u)) >= MACHEPS*zv_norm2(x)*cond_est)
  479.     {
  480.     errmesg("zLUfactor()/zLUsolve()");
  481.     printf("# LU solution error = %g [cf MACHEPS = %g]\n",
  482.            zv_norm2(zv_sub(x,z,u)), MACHEPS);
  483.     }
  484.  
  485.     zvm_mlt(B,z,y);
  486.     zv_copy(y,x);
  487.     tracecatch(zLUAsolve(A,pivot,x,x),"main");
  488.     if ( zv_norm2(zv_sub(x,z,u)) >= MACHEPS*zv_norm2(x)*cond_est)
  489.     {
  490.     errmesg("zLUfactor()/zLUAsolve()");
  491.     printf("# LU solution error = %g [cf MACHEPS = %g]\n",
  492.            zv_norm2(zv_sub(x,z,u)), MACHEPS);
  493.     }
  494.  
  495.     MEMCHK();
  496.  
  497.     /* QR factorisation */
  498.     zm_copy(B,A);
  499.     zmv_mlt(B,z,y);
  500.     notice("QR factor/solve:");
  501.     diag = zv_get(A->m);
  502.     zQRfactor(A,diag);
  503.     zQRsolve(A,diag,y,x);
  504.     if ( zv_norm2(zv_sub(x,z,u)) >= MACHEPS*zv_norm2(x)*cond_est )
  505.     {
  506.     errmesg("zQRfactor()/zQRsolve()");
  507.     printf("# QR solution error = %g [cf MACHEPS = %g]\n",
  508.            zv_norm2(zv_sub(x,z,u)), MACHEPS);
  509.     }
  510.     printf("# QR cond(A) approx= %g\n", zQRcondest(A));
  511.     Q = zm_get(A->m,A->m);
  512.     zmakeQ(A,diag,Q);
  513.     zmakeR(A,A);
  514.     zm_mlt(Q,A,C);
  515.     zm_sub(B,C,C);
  516.     if ( zm_norm1(C) >= MACHEPS*zm_norm1(Q)*zm_norm1(B) )
  517.     {
  518.     errmesg("zQRfactor()/zmakeQ()/zmakeR()");
  519.     printf("# QR reconstruction error = %g [cf MACHEPS = %g]\n",
  520.            zm_norm1(C), MACHEPS);
  521.     }
  522.  
  523.     MEMCHK();
  524.  
  525.     /* now try with a non-square matrix */
  526.     A = zm_resize(A,15,7);
  527.     zm_rand(A);
  528.     B = zm_copy(A,B);
  529.     diag = zv_resize(diag,A->n);
  530.     x = zv_resize(x,A->n);
  531.     y = zv_resize(y,A->m);
  532.     zv_rand(y);
  533.     zQRfactor(A,diag);
  534.     x = zQRsolve(A,diag,y,x);
  535.     /* z is the residual vector */
  536.     zmv_mlt(B,x,z);    zv_sub(z,y,z);
  537.     /* check B*.z = 0 */
  538.     zvm_mlt(B,z,u);
  539.     if ( zv_norm2(u) >= 100*MACHEPS*zm_norm1(B)*zv_norm2(y) )
  540.     {
  541.     errmesg("zQRfactor()/zQRsolve()");
  542.     printf("# QR solution error = %g [cf MACHEPS = %g]\n",
  543.            zv_norm2(u), MACHEPS);
  544.     }
  545.     Q = zm_resize(Q,A->m,A->m);
  546.     zmakeQ(A,diag,Q);
  547.     zmakeR(A,A);
  548.     zm_mlt(Q,A,C);
  549.     zm_sub(B,C,C);
  550.     if ( zm_norm1(C) >= MACHEPS*zm_norm1(Q)*zm_norm1(B) )
  551.     {
  552.     errmesg("zQRfactor()/zmakeQ()/zmakeR()");
  553.     printf("# QR reconstruction error = %g [cf MACHEPS = %g]\n",
  554.            zm_norm1(C), MACHEPS);
  555.     }
  556.     D = zm_get(A->m,Q->m);
  557.     zmam_mlt(Q,Q,D);
  558.     for ( i = 0; i < D->m; i++ )
  559.     m_set_val(D,i,i,zsub(m_entry(D,i,i),ONE));
  560.     if ( zm_norm1(D) >= MACHEPS*zm_norm1(Q)*zm_norm_inf(Q) )
  561.     {
  562.     errmesg("QRfactor()/makeQ()/makeR()");
  563.     printf("# QR orthogonality error = %g [cf MACHEPS = %g]\n",
  564.            zm_norm1(D), MACHEPS);
  565.     }
  566.  
  567.     MEMCHK();
  568.  
  569.     /* QRCP factorisation */
  570.     zm_copy(B,A);
  571.     notice("QR factor/solve with column pivoting");
  572.     pivot = px_resize(pivot,A->n);
  573.     zQRCPfactor(A,diag,pivot);
  574.     z = zv_resize(z,A->n);
  575.     zQRCPsolve(A,diag,pivot,y,z);
  576.     /* pxinv_zvec(pivot,z,x); */
  577.     /* now compute residual (z) vector */
  578.     zmv_mlt(B,x,z);    zv_sub(z,y,z);
  579.     /* check B^T.z = 0 */
  580.     zvm_mlt(B,z,u);
  581.     if ( zv_norm2(u) >= MACHEPS*zm_norm1(B)*zv_norm2(y) )
  582.     {
  583.     errmesg("QRCPfactor()/QRsolve()");
  584.     printf("# QR solution error = %g [cf MACHEPS = %g]\n",
  585.            zv_norm2(u), MACHEPS);
  586.     }
  587.  
  588.     Q = zm_resize(Q,A->m,A->m);
  589.     zmakeQ(A,diag,Q);
  590.     zmakeR(A,A);
  591.     zm_mlt(Q,A,C);
  592.     ZM_FREE(D);
  593.     D = zm_get(B->m,B->n);
  594.     /******************************
  595.     px_cols(pivot,C,D);
  596.     zm_sub(B,D,D);
  597.     if ( zm_norm1(D) >= MACHEPS*zm_norm1(Q)*zm_norm1(B) )
  598.     {
  599.     errmesg("QRCPfactor()/makeQ()/makeR()");
  600.     printf("# QR reconstruction error = %g [cf MACHEPS = %g]\n",
  601.            zm_norm1(D), MACHEPS);
  602.     }
  603.     ******************************/
  604.  
  605.     /* Now check eigenvalue/SVD routines */
  606.     notice("complex Schur routines");
  607.     A = zm_resize(A,11,11);
  608.     B = zm_resize(B,A->m,A->n);
  609.     C = zm_resize(C,A->m,A->n);
  610.     D = zm_resize(D,A->m,A->n);
  611.     Q = zm_resize(Q,A->m,A->n);
  612.  
  613.     MEMCHK();
  614.  
  615.     /* now test complex Schur decomposition */
  616.     /* zm_copy(A,B); */
  617.     ZM_FREE(A);
  618.     A = zm_get(11,11);
  619.     zm_rand(A);
  620.     B = zm_copy(A,B);
  621.     MEMCHK();
  622.  
  623.     B = zschur(B,Q);
  624.     MEMCHK();
  625.  
  626.     zm_mlt(Q,B,C);
  627.     zmma_mlt(C,Q,D);
  628.     MEMCHK();
  629.     zm_sub(A,D,D);
  630.     if ( zm_norm1(D) >= MACHEPS*zm_norm1(Q)*zm_norm_inf(Q)*zm_norm1(B)*5 )
  631.     {
  632.     errmesg("zschur()");
  633.     printf("# Schur reconstruction error = %g [cf MACHEPS = %g]\n",
  634.            zm_norm1(D), MACHEPS);
  635.     }
  636.  
  637.     /* orthogonality check */
  638.     zmma_mlt(Q,Q,D);
  639.     for ( i = 0; i < D->m; i++ )
  640.     m_set_val(D,i,i,zsub(m_entry(D,i,i),ONE));
  641.     if ( zm_norm1(D) >= MACHEPS*zm_norm1(Q)*zm_norm_inf(Q)*10 )
  642.     {
  643.     errmesg("zschur()");
  644.     printf("# Schur orthogonality error = %g [cf MACHEPS = %g]\n",
  645.            zm_norm1(D), MACHEPS);
  646.     }
  647.  
  648.     MEMCHK();
  649.  
  650.     /* now test SVD */
  651.     /******************************
  652.     A = zm_resize(A,11,7);
  653.     zm_rand(A);
  654.     U = zm_get(A->n,A->n);
  655.     Q = zm_resize(Q,A->m,A->m);
  656.     u = zv_resize(u,max(A->m,A->n));
  657.     svd(A,Q,U,u);
  658.     ******************************/
  659.     /* check reconstruction of A */
  660.     /******************************
  661.     D = zm_resize(D,A->m,A->n);
  662.     C = zm_resize(C,A->m,A->n);
  663.     zm_zero(D);
  664.     for ( i = 0; i < min(A->m,A->n); i++ )
  665.     zm_set_val(D,i,i,v_entry(u,i));
  666.     zmam_mlt(Q,D,C);
  667.     zm_mlt(C,U,D);
  668.     zm_sub(A,D,D);
  669.     if ( zm_norm1(D) >= MACHEPS*zm_norm1(U)*zm_norm_inf(Q)*zm_norm1(A) )
  670.     {
  671.     errmesg("svd()");
  672.     printf("# SVD reconstruction error = %g [cf MACHEPS = %g]\n",
  673.            zm_norm1(D), MACHEPS);
  674.     }
  675.     ******************************/
  676.     /* check orthogonality of Q and U */
  677.     /******************************
  678.     D = zm_resize(D,Q->n,Q->n);
  679.     zmam_mlt(Q,Q,D);
  680.     for ( i = 0; i < D->m; i++ )
  681.     m_set_val(D,i,i,m_entry(D,i,i)-1.0);
  682.     if ( zm_norm1(D) >= MACHEPS*zm_norm1(Q)*zm_norm_inf(Q)*5 )
  683.     {
  684.     errmesg("svd()");
  685.     printf("# SVD orthognality error (Q) = %g [cf MACHEPS = %g\n",
  686.            zm_norm1(D), MACHEPS);
  687.     }
  688.     D = zm_resize(D,U->n,U->n);
  689.     zmam_mlt(U,U,D);
  690.     for ( i = 0; i < D->m; i++ )
  691.     m_set_val(D,i,i,m_entry(D,i,i)-1.0);
  692.     if ( zm_norm1(D) >= MACHEPS*zm_norm1(U)*zm_norm_inf(U)*5 )
  693.     {
  694.     errmesg("svd()");
  695.     printf("# SVD orthognality error (U) = %g [cf MACHEPS = %g\n",
  696.            zm_norm1(D), MACHEPS);
  697.     }
  698.     for ( i = 0; i < u->dim; i++ )
  699.     if ( v_entry(u,i) < 0 || (i < u->dim-1 &&
  700.                   v_entry(u,i+1) > v_entry(u,i)) )
  701.         break;
  702.     if ( i < u->dim )
  703.     {
  704.     errmesg("svd()");
  705.     printf("# SVD sorting error\n");
  706.     }
  707.     ******************************/
  708.  
  709.     ZV_FREE(x);    ZV_FREE(y);    ZV_FREE(z);
  710.     ZV_FREE(u);    ZV_FREE(diag);
  711.     PX_FREE(pi1);    PX_FREE(pi2);    PX_FREE(pivot);
  712.     ZM_FREE(A);    ZM_FREE(B);    ZM_FREE(C);
  713.     ZM_FREE(D);    ZM_FREE(Q);
  714.  
  715.     mem_stat_free(1);
  716.  
  717.     MEMCHK();
  718.     printf("# Finished torture test for complex numbers/vectors/matrices\n");
  719.     mem_info();
  720. }
  721.  
  722.